package com.zznode.client.wlanpolling.utils;

public class GaussArithmetic {
	private static final double EPSILON = 1e-10;
	
	/**
	 * for test.
	 * 
	 * @param args
	 */
	public static void main(String[] args) {
		double[][] a = { { 2, 0, 0 }, { 3, 1, 0 }, { 0, 4, 1 } };
		double[] b = { 2, 4, 5 };
		double[] x = { 0, 0, 0 };
		int n;
		n = b.length;

		if (GaussEliminate(a, b, x)) {
			for (int i = 0; i < n; i++) {
				System.out.println(x[i]);
			}
		} else {
			System.out.println("This is a single matrix!");
		}
	}

	public static boolean GaussEliminate(double[][] a, double[] b, double[] x) {
		// Gaussian elimination method
		int maxI, i, j, n, n1;
		double maxV, s1;

		n = b.length;
		n1 = n - 1;
		
		for (int ii = 0; ii < n; ii++) {
			maxV = Math.abs(a[ii][ii]);
			maxI = ii;
			for (i = ii + 1; i < n; i++) {
				if (Math.abs(a[i][ii]) > maxV) {
					maxV = Math.abs(a[i][ii]);
					maxI = i;
				}
			}

			if (maxV == 0.0) {
				return (false);
			}

			if (ii != maxI) {
				for (j = ii; j < n; j++) // exchange rows
				{
					s1 = a[maxI][j];
					a[maxI][j] = a[ii][j];
					a[ii][j] = s1;
				}
				s1 = b[maxI];
				b[maxI] = b[ii];
				b[ii] = s1;
			}
			
			for (i = ii + 1; i < n; i++) // elemination
			{
				s1 = a[i][ii] / a[ii][ii];
				for (j = ii + 1; j < n; j++) {
					a[i][j] -= a[ii][j] * s1;
				}
				b[i] -= b[ii] * s1;
			}
		}
		
		x[n1] = b[n1] / a[n1][n1]; // back substitution
		for (i = n1; i >= 0; i--) {
			s1 = 0.;
			for (j = i + 1; j < n; j++) {
				s1 += a[i][j] * x[j];
			}
			x[i] = (b[i] - s1) / a[i][i];
		}
		
		return true;
	}
	
    public static double[] lsolve(double[][] A, double[] b) {
        int N  = b.length;

        for (int p = 0; p < N; p++) {

            // find pivot row and swap
            int max = p;
            for (int i = p + 1; i < N; i++) {
                if (Math.abs(A[i][p]) > Math.abs(A[max][p])) {
                    max = i;
                }
            }
            double[] temp = A[p]; A[p] = A[max]; A[max] = temp;
            double   t    = b[p]; b[p] = b[max]; b[max] = t;

            // singular or nearly singular
            if (Math.abs(A[p][p]) <= EPSILON) {
                throw new RuntimeException("Matrix is singular or nearly singular");
            }

            // pivot within A and b
            for (int i = p + 1; i < N; i++) {
                double alpha = A[i][p] / A[p][p];
                b[i] -= alpha * b[p];
                for (int j = p; j < N; j++) {
                    A[i][j] -= alpha * A[p][j];
                }
            }
        }

        // back substitution
        double[] x = new double[N];
        for (int i = N - 1; i >= 0; i--) {
            double sum = 0.0;
            for (int j = i + 1; j < N; j++) {
                sum += A[i][j] * x[j];
            }
            x[i] = (b[i] - sum) / A[i][i];
        }
        return x;
    }

}
